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Abstract 

We examine the approach to equilibrium of the micromaser. Analytic methods 
are first used to show that for large times (i.e. many atoms) the convergence is 
governed by the next to leading eigenvalue of the corresponding discrete evolution 
matrix. The model is then studied numerically. The numerical results confirm 
the phase structure expected from analytic approximation methods and agree for 
large times with the analysis of Elmfors et al in terms of the "continuous master 
equation". For short times, however, we see evidence for interesting new structure 
not previously reported in the literature. 
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1 Introduction 



The micromaser|l| provides an excellent theoretical and experimental testing 
ground for many fundamental properties of cavity quantum electrodynamics 
and quantum mechanics in general. The physical situation under consider- 
ation consists of a superconducting, high Q cavity, that is being traversed 
by a low intensity beam of two state atoms. The atoms interact with the 
electromagnetic field in the cavity via an electric dipole interaction. The dy- 
namics of the atom cavity system is well described by the Jaynes-Cummings 
model0 . If the cavity transit time r is short compared to the average time 
T between atoms, there is effectively only one atom in the cavity at a time 
and the atoms in the beam interact with each other only via their residual 
effect on the electromagnetic field. For example if the atoms enter the cavity 
preferentially in their excited state, and emit a photon, the photons tend 
to build up inside the cavity, and each successive atom sees a stronger pho- 
ton field when it enters the cavity. This "pumping" is responsible for the 
evolution of the system into a microwave laser, or "maser" . 

Three independent time scales determine the overall dynamical behaviour: 
the time interval T between consecutive atoms, r the time spent by each atom 
inside the cavity and I/7, the characteristic photon decay time 7 inside the 
cavity. An important physical quantity is the dimensionless "pumping rate" 
N = R/j, where R = 1/T and 7 is the characteristic photon decay time of 
the cavity. N can be thought of as the number of atoms that pass through the 
cavity in a single photon decay time. When both damping and pumping are 
present, the photon distribution inside the cavity asymptotically approaches 
a steady state distribution. The details of the steady state (equilibrium) 
distribution depend on the time r that the atom spends in the cavity as well 
as the dimensionless pumping rate N. 

Although much work has been done on the equilibrium properties of this 
system, to the best of our knowledge there has not been a systematic anal- 
ysis of the initial stages of the approach to equilibrium, which in principle 
can be important in determining the outcome of very low flux experiments. 
The purpose of the present work is to examine numerically this approach to 
equilibrium. In particular, we will see how varying the physical parameters 
affects the rate at which equilibrium is reached: i.e. how many atoms must 
pass through the cavity before a steady state photon distribution is estab- 
lished. In a recent paper, Elmfors et aZ0 looked at long time correlations 
in the outgoing atomic beam and their relation to the various phases of the 
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micromaser system. The properties they considered were associated with 
the equilibrium configuration of the cavity photon distribution, but there is 
a close connection between the correlation functions considered in [[J and 
the near equilibrium dynamical behaviour that we will be examining. As we 
will show, our results agree with those of [f| in the appropriate (i.e large N) 
limits. 

The paper is organized as follows: In Section 2 we review the JC model 
and its application to the physical situation at hand. In particular we de- 
rive transition matrix S that governs the master equation for the dynamical 
evolution of the photon distribution inside the cavity. We also derive the 
expression for the probability P(+) of finding the atom in the excited state. 
In Section 3 we show that the approach to equilibrium of the photon dis- 
tribution and of the physically measured P(+) is governed by the leading 
eigenvalues of S. We compare our theoretical analysis to that of Elmfors 
et al, who looked at correlation functions instead of the photon distribution 
directly. In Section 4, we describe the numerical experiment that we use to 
analyze the approach to equilibrium, and compare our results to our theoret- 
ical analysis and to that of Elmfors et al. Section 5 closes with a summary 
and conclusions. 



2 The Jaynes-Cummings Model 

We consider atoms with two possible states |±) with energy difference 

E+ -E^ = hu a (1) 

For a high Q cavity, the electromagnetic field is well approximated by a single 
mode, with energy E c = hu> c . For simplicity we assume that the cavity is 
tuned so that its fundamental frequency is equal to that of the atom: 

U a = L0 C = LO (2) 

For a single atom traversing the cavity, the dynamics of the atom-cavity 
system is governed by the JC Hamiltonian. 

H = ua^a H — uoa z + g{aa + + a)a_) (3) 

where g is the coupling constant, (a) are the photon creation (annihilation) 
operators and a± = i^B^v) are operators which raise and lower the atomic 
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states (cr x , <j y , and a z are the Pauli matrices). In the absence of the dipole 
interaction (i.e. when g = 0) the atom-plus-field energy eigenstates are \n, s), 
where n = 0, 1, ... is the the photon number and s — ± for the two atomic 
levels. When g is non-zero the system makes transitions between the energy 
eigenstates of the non-interacting system with probabilities, 

|(n,-|e-^|n,-)| 2 = \-q n {r) 

|(n-l,+|e-^K-)| 2 = q n (r) (4) 

\(n,+\e- mt \n,+)\ 2 = 1 - q n+1 (r) 

|(n + l,-|e-^K+)| 2 = q n+1 (r) 

These probabilities are expressed in terms of the quantity 

q n (r) = sin 2 (g^/nrj (5) 

This is a completely solvable quantum mechanical system. We suppose that 
the atom/ cavity states are uncorrelated at t — 0, so that 

= \lpatom) ® \lpcav) 

= («|+>+)9|-»®(EC»|n» ( 6 ) 

n 

The interaction between the atom and electromagnetic field causes the states 
to be entangled. The exact result for the wave function after an interaction 
time t is: 

1^(0) = [( aC n cos(gVn + It) - if3C n+1 sin(gVn + It)) \n, +) 

n 

+ i^—iaCn^i s\n{gy/nt) + j3C n cos (g^/nt)j \n, — ) (7) 

We now define P n , s (t) as the probability of finding the atom in the state s, 
and n photons in the cavity. Specifically one has: 

P n , + (t) = (n,+|^(t))| 2 

= aP n (l- q n+ i(t)) + bP n+1 q n+1 (t) 

— [ a (l ~~ Qn+l)^n,m + ^n+l^m,n+l]^m 

=: M(t,+) nmPm (8) 

P n ,_(t) = |(n,-|^(t))| 2 

= aP n _ iqn (t) + bP n (l- q n (t)) 

= [aq n S m , n -i + 6(1 ~ <?n)<5„, m ]P m 

=: M(t,-) nm P m (9) 
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where a = a*a (b = (3*(3) is the probability that the atom entered the 
cavity in the excited (lower) state, while P n = C*C n is the probability that 
there were n photons in the cavity initially It follows directly that the 
probabilities V+(t) and V-{t) of finding the atom in the upper and lower 
states, respectively, for unknown cavity state, are: 

r + (t) = £K^+hH*)>l 2 

n 

= J2( aP n(~L-<ln+l(t))+bP n+1 q n+1 (t)) (10) 

n 

v-{t) = £|<n,-|V>(*)>l 2 

n 

= J2( aP n-Mt) + bP n (l-q n (t))) (11) 

n 

Eqs.([ll|) and ([HI can be written in matrix form: 

V s (t) = J2M(t,s) nm P m (12) 

n,m 

where M(t,s) nm is defined in Eqs(Eq.(||)) and (Eq.(|)) above. 

Conversely, if we are not interested in determining the state of the atom 
then the probability of finding exactly n photons in the cavity isQ: 

Pn{t) = P n ,+ {t) + Pn,-{t) 

= M(t) nm P m (13) 



where 



+ bq n+1 + (1 - aq n+ i - bq n )8 n)m (14) 



Eqs.([T3|) and ( |I4]) give the master equation for the time evolution of the 
photon distribution in the presence of the atom-cavity interaction, without 
thermal dissipation. The first term in the transition matrix M gives the 
probability that an n-photon state occurs through decay of the excited atomic 
state in interaction with a cavity containing n — 1 photons. It is given by the 
product of a (the probability for the atom to be in the excited state) times 

: We henceforth adopt the convention that repeated indices are to be summed, unless 
stated otherwise. 
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p n -i (the probability that the cavity contains n — 1 photons) times q n (r) 
(probability for a transition between an unperturbed eigenstate \n — 1, +) 
and an unperturbed state \n, — ). Similarly, the second term comes from the 
excitation of the atomic ground state, and the third term is the contribution 
from processes that leave the atom unchanged. 

The above analysis assumes that the system under consideration is in 
a pure quantum state. In a realistic experiment both the atom and cavity 
would be described by a density matrix representing a mixed state. We will 
now show that under some simple assumptions the above formulas also apply 
to the more realistic case. Let p{t) a c denote the density matrix describing 
the atom/cavity system at time t. Operator expectation values are given by, 

(6) = Tr aC (p(t)acO) (15) 

where the subscript aC indicates that the trace is over both the atomic and 
the cavity states. If we restrict ourselves to the measurement of observables 
that involve only atomic operators, the expectation value of such an operator 
is given by, 

(6(t)) = Tr a (d(t)Tr c p(t)ac) = Tr a (6(t)p(t)red) (16) 
where the operator 

p{t) red = Tr c p{t) (17) 

is the trace over the cavity states of the total density matrix and is called the 
reduced density matrix of the system. To determine the expectation values 
of atomic operators at arbitrary times we need to know p(t) re d for any time 
t. 

Our system consists of a series of atoms that enter a cavity containing 
electromagnetic radiation. We require that the time T between atoms, and 
the photon decay time I/7, are much larger than the time that any given 
atom spends in the cavity (T ^> r, I/7 ^> r); equivalently we assume that 
the time scale of interactions within the reservoir is much smaller than the 
time scale over which we want to consider the evolution of the system. Under 
this condition the density matrix for the initial state of the system can be 
factored into a product of density matrices for the cavity and the individual 
atoms: 

P = PC®Pa (18) 
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Substituting into Eq.(|H]) we have, 



p(t) red = Tr c (p(t) a <8> p(t) c ) (19) 

We treat the cavity as a reservoir and sum over the large number of reservoir 
states to obtain the ensemble averaged density matrix, 

■y M oo 

Pc ^ Pc = ^M^ooj^^Zpn = ^Pn\n){n\ (20) 



n=l 



where p n > and S^L p n — 1. We assume that the incoming atoms are 
uncorrelated and have initial states that can be represented as a diagonal 
mixture of excited and unexcited states with a density matrix of the form, 

A>=(q °)=a|+><+| + &|->H (21) 

where a, b > and a + b = 1. Physically, a is the probability that the 
atom is initially in the excited state, and b is the probability that the atom 
is initially in the ground state. Combining Eq.(ppp and Eq.(^Tj) we have 
an expression for the atom-cavity density matrix at the time the atom first 
enters the cavity: 

p(0) = J2(aPn\n+)(n+\ + bp n \n-)(n-\) (22) 

n 

To study the time evolution of the atomic variables, we need the time depen- 
dent reduced density matrix at time t. A straightforward calculation reveals 
that 



p(t) red = Tr c (e- im p(0)e 

( r+(t) o 

V V-{t) 



iHt 



(23) 



with V s (t) given in Eq. ([T0|) and Eq . (|TT]) . Similarly, if we are interested in 
measuring only cavity observables, we must consider the reduced density 
matrix for the cavity, which after time, t is given by: 



JHt 



PcMt) = Tr a (e^'p(0)e l 

= J2 p n(t)\n)(n\ (24) 
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where P n (t) is given by Eq.(|T^). Thus if the initial atomic and cavity density 
matrices are diagonal, then they both remain diagonal, and they give rise to 
precisely the same master equation for the photon probability distribution 
as in the case of pure states. 

Eq.flLll) can be modified to include the effects of thermal damping. Sup- 
pose the photon distribution inside the cavity initially is p(°\ An atom enters 
the cavity and exits after an interaction time r. Assuming that 7T <C 1, we 
neglect damping during the atom-photon interaction. The probability distri- 
bution for the photons just before the next atom enters the cavity is@: 

p(T) = e - 7LcT M(r)p (0) (25) 

where T is the time between atoms and 



{Lc)nm = {nb + l)(n5 n , m — (n + l)<5 n+ i jm ) + n b ((n + 1)5 

n,m Tl&n— l,m ) (26) 

We can also take into account the fact that atoms in the beam arrive at time 
intervals that are Poisson distributed, with an average time interval of T = 
1/R between them. Multiplying by the distribution function exp(— RT) RT 
and integrating we find the averaged photon distribution just prior to the 
arrival of the second atom to be: 

(p(T)) T = Sp<® (27) 

S = — : — M (28) 

l + L c /N V ; 

This is the form of the master equation that we will use to describe the 
dynamics of the photon distribution inside the cavity. We will refer to it as 
the discrete master equation. 

The analysis in || starts from a different master equation, which is called 
the continuous master equation. To derive the continuous master equation, 
Elmfors et al consider a situation where the flux of incoming atoms is large 
enough that the atoms have Poisson distributed arrival times, so that each 
atom has the same probability Rdt of arriving in an infinitesimal time dt. 
They further assume that the interaction within the cavity takes place in 
a time much less than this time interval (r dt) which means that the 
interaction is essentially instantaneous. The contributions from damping 
and pumping during the time interval d t can then be considered separately. 
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The contribution from the damping is exactly the same as before [Eq.(^])]. 
The contribution from pumping has the form, 

{dp)„={M-l)Rpdt (29) 

where M is given in Eq . (|14|) as before. The continuous form of the master 
equation is obtained by combining the two contributions, 

dp = {—^L c p + (M — l)Rp}dt =: —'jLpdt (30) 
L = (n b + 1) {n5 n;m - (n+ 1) 5„+i, m ) + n fe ((™ + 1) 5„, m - n5„, m+ i) 
+ AT ((og n+1 + 6g„) 5 n , m - ag n 5 n>m+ i - &g n +i<Wi,m) (31) 

It is expected that the continuous and discrete master equations agree when 
the number of atoms per photon decay time is very large.0. In particular, 
when the thermalization time scale I/7 is much greater than the time T 
between atoms, the large time (many atom) dynamics should agree. This 
correspondence will be verified explicitly below. 



3 Eigenvalue Problem and Approach to Equi- 
librium 

Recall the physical system we are considering. We inject a series of atoms 
into a cavity. The time, T, between the atoms is much greater than the 
time an individual atom takes to pass through the cavity r so that there 
is never more than one atom in the cavity at one time. Each atom inter- 
acts with the cavity's electromagnetic field as it passes through the cavity, 
and consequently, the photon distribution changes repeatedly. According to 
the discrete master equation Eq.(|27|) after k atoms have passed through the 
cavity, the photon distribution, p(k) is: 

p(k) = S k p(0) } (32) 

where p(0) is the initial photon distribution. Equilibrium occurs when the 
photon distribution is no longer changed by the transition matrix S. That is 

Sp eq = p eq (33) 

2 Note that the atomic flux must still be small enough so that effectively only one atom 
is in the cavity at a time. 
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and the equilibrium photon distribution p eq is an eigenvector of S with eigen- 
value 1. We expect that as k — > oo, p(k) — > p eq , and it is precisely this 
approach to equilibrium that we wish to investigate. 

As discussed in |3j the approach to equilibrium is governed by the eigen- 
values of S. Before proving this, we need some preliminaries. The right 
eigenvectors of the matrix S [Eq.(|28|)l are written pW, the left eigenvectors 
are u T ^ and the eigenvalues are K^ l \ 

Sp® = «V (34) 

u T(l) S = u T(l) K (l) (35) 

Eq.(|33"D implies that there is an eigenvector with eigenvalue unity. As we will 
see, all other eigenvalues must be less than one in order for the system to be 
stable. For convenience, we label the eigenvalues by size; 

«W = 1 > > k®.... (36) 

With this labelling, p^ 1 ' = p eq and 

Sp w = p {1) (37) 
U T(D S = u T(i) (38) 

From Eq.(Trj), Eq.(|2B[) and Eq.(B5[) it follows that u T ^ is a vector with all 
components equal to one. To see this, note that the transformation matrix S 
must preserve the norm of the probability distribution. It therefore follows 
that 

Yl SnmPm = ^Pn (39) 
n,m n 

for all p n , which in turn requires: 

J2 S nm=l (40) 

n 

for all m. This is equivalent to Eq.(|38|) with u^ 1 ^ = 1. 

We will use the fact that the similarity transform of the form 

Tml = pS (41) 

diagonalizes the matrix S and that the inverse of T is given by the left 
eigenvector: 

Tim = u T ® (42) 



9 



Finally, from Eq.([IT]) and Eq.([K|) we have, 

U m l Pm = TimT ma = ^la (43) 

These results allow us to write, the evolution matrix as, 

S nm = (44) 

It is easy to verify that this expression satisfies the eigenvector equations 
Eq.(|3§ and Eq.(£§). 

In order to investigate the approach to equilibrium we start with the fact 
that 

p(k + 1) = Sp(k) (45) 

Now define 

dp{k) = p{k + 1) - p{k) = {S- l)p{k)dk (46) 

where dk = {k + 1) — k — 1. This is the discrete version of the continuous 
evolution equation Eq. (|30|) . In the limit of large k, it can be treated as a 
differential equation. In order to integrate it, we first define 

Q(k) = T-'pik) (47) 

so that Eq.(^) reads in component form: 

dQn(k) = (/t (n) - l)Q n (k)dk (48) 

Eq.(f48|) can be trivially integrated to give 

Qn(k) = Q n (0)e-^^ k (49) 
The solution for the photon distribution after the kth atom is therefore: 
p n (k) = ^T nm exp(-(l- K ( m ))A:)T m ^(0) 

m,l 

= Y.PneM-a-K {m) )k)u^ m) Pl (0) (50) 

m,l 

where we have used Eq . (|4l|) and Eq . (^2|) . As k — > oo only the leading 
eigenvalue = 1 survives, and determines the asymptotic value of p n (k). 
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The next to leading eigenvalue will control the rate of convergence. In 
particular: 



p n (k) ^p n (oo)+^exp(-(l-KW)A:)(Ap)« (51) 

l>2 

where we have defined 



(Ap)P=pP£^W0) (52) 

m 

and 

Pn(00) = P«E«r (1) Pl(0) 
/ 

= P ( n ] (53) 

Thus, the photon distribution converges to its equilibrium value, as desired. 
Moreover, Eq.([H]) implies: 

p n (k) - Pn (oo) = £exp(-(l - ^)k)(Ap)V) » exp(-(l - ^)(Ap)( 2 ) + . . 

l>2 

and the approach to equilibrium is determined by if there is no degen- 
eracy in the next to leading eigenvalues. As we will see in the subsequent 
numerical analysis, interesting things occur when and are close within 
the numerical accuracy of the calculation. 

Note that the continuous master equation Eq.(^) can be integrated in 
precisely the same way, with —(1 — S) replaced by — jL, so that the asymp- 
totic behaviour is controlled by the eigenvalues of — jL, instead of 1 — k®. 
In the appendix we prove that these eigenvalues coincide in the large N limit, 
and this is verified numerically in the next section. 

Before closing this section we make one further comparison between our 
analysis and that of 0, who look at correlation functions of the spin variables 
of the form 

(ss)k- {s} 2 

^ k = Tlx Tf 55 

(s 2 ) - {s 2 



where 



( S ) = £ sV(s); V{s) = £ S(s) nmP % = ul^S{s) nm p^ (56) 

s=±l n 
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Similarly, (ss)k is the joint probability for observing the states of two atoms, 
Si followed by s 2 , with k unobserved atoms between them: 



(ss) k 



s 1 s 2 V k (s 1 .s 2 ) 



S1=±1,S2=±1 



Vk(si,s 2 ) 



u T ^S(s 2 )S k S( Sl )p^ 



(57) 



Note that the above expectation values have assumed that the photon distri- 
bution is already at its equilibrium value before the spin of the first atom 
is measured. In spite of this, the correlation function Eq.flSop does describe 
the approach to equilibrium and is directly related to the quantities that we 
study in the present paper. In particular, the correlation length associated 
with 7^ approaches zero in the limit of large k, and this approach is gov- 
erned by the same eigenvalues as the approach to equilibrium of the photon 
distribution. The basic argument is as follows. In effect (ss)fc depends on 
the conditional probability that one first measures the spin to be s\, say and 
that after k atoms pass, spin s 2 is measured. However, when one measures 
the first spin, one effectively applies a projection operator which moves the 
photon distribution away from the equilibrium configuration. The shape of 
this projected photon distribution then determines the correlation between 
the first and second spin measurements. One expects that for a very large 
number k of atoms between the two measurements, the photon distribution 
again approaches its equilibrium value, so that correlation between the two 
spin measurement vanishes. That is: 



and the correlation function approaches zero as k approaches infinity. The 
correlations between well separated atoms therefore measure the rate at 
which the photon distribution settles back to equilibrium. 

We will now prove the above assertions. We look for exponential decay 
of the correlation function, 



where k ~ Rt and ( is the correlation length, or the typical length of time 
over which the cavity remembers pumping events. We can rewrite, 



lim fc _ 00 (ss) fc -> (s) 2 



(58) 




(59) 



1 



TD k T~ 



1 



(60) 
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where 

D kl = K k 5 kl . (61) 

We want to consider the way in which D k approaches equilibrium. We start 
from, 

D k+i _ D k = p _ ^ D k ^ 

We consider the limit of large k in order to isolate the exponential depen- 
dence. We write, D k+1 - D k -> dD(k) and 1 = (k + 1) - k -> dk. We 
obtain, 

dD nm (k) = (D - l) nl D lm (k) dk 

= (k^ -l)5 nl D lm (k) dk (63) 
= (k^ -l)D nm (k) dk 

which has the solution, 

D nm {k) = D nm (0)exp-( 1 - K(n) ) fc (64) 

where D nm (0) is an integration constant. Writing out the first two terms in 
the sum over n we obtain, 

D nm {k) = 5 nl D nrn (0)ex p -^-^ k + S^ =2J D nm (0)exp-( 1 - K( " ) ) fe (65) 

Since = 1 and D is diagonal we have, 

D nm (k) = 5 nl 5 ml D nm (0) + S^ =2 D nm (0)exp-( 1 - K(n) ) fc (66) 

Substituting Eq . (|60D and Eq . (|66"D into Eq.(^7]) we have, 

(ss) k = S S1iS2 sis 2 m^ (1) S nj (s 2 )T jq 

S ql 5 ml D qm (0) + S^ o =2J D (2m (0)exp-( 1 - K(9) ) fe ] T-lSn^p? 

Denoting the first (leading order) term in the square brackets by (ss) z fc '°', we 
find 

(88)t = S Sl , S2SlS2 ^ 1 )^( S2 )T ilJ D 11 (0)Tr r 1 Sw( S i)ri 1) 

= E Sl>S2 s 1 s 2 u^ 1) S nj (s 2 )pf ) D 11 (0)u^S rl (s 1 )p < i 1) (67) 



13 



where we have used Eq.(¥T) and Eq.(f4"2|). Since /c' 1 ' = 1 we have from Eq.| 
that Du is independent of k and from Eq.fl6T|) that Dn = 1. This gives, 



The second term in square brackets of Eq.(|67|) shows that, as stated earlier, 
all of the eigenvalues other than must be less than one, or the correla- 
tion length diverges. Recalling that the eigenvalues are labelled by size, the 
leading order non-zero term in the correlation function Eq. (|55|) has the form, 



Ik ~ e~ (1 ~ K(2))fc (69) 



and from Eq.(|59|) we have, 



Comparing Eq.(|51~D and Eq.(^) we see that the correlation length ( is deter- 
mined by the same eigenvalue that determines the approach to equilibrium, 
as claimed. 



4 Numerical Results 

We wish to investigate numerically the approach to equilibrium of the photon 
distribution as described by the dynamical master equation Eq.(^), with S 
given by Eq.(^). We assume for concreteness that before the first atom 
enters the cavity, the photons are in thermal equilibrium at the tempera- 
ture T which characterizes the thermal properties of the system throughout 
the experiment. In principle the asymptotic properties of the approach to 
equilibrium should not be sensitive to the initial photon distribution, but 
we observed some interesting short time behaviour which presumably does 
depend on the initial distribution. The short time behaviour may thus have 
physical relevance. To begin, we start with the photon distribution: 

p n = [1 - e kT] e kT (71) 

The mean photon number is therefore 

n b = Y, n np n = (72) 

ei=T — 1 
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Note that the mean photon number also plays an important role in the master 
equation, i.e. in the matrix Lq (Eq.fl26|)) that determines the rate at which 
the photon distribution relaxes into thermal equilibrium. For a typical two 
level Rydberg atom ||, hoo « 1.4 x 10~ 23 J, so that 

T(Kelvin) « 1 —— 

v ; In 1 + l/n b 

Typical experimental temperatures range between T = OAK (rib = 0.1) and 
T = 10K (rib = 10). In terms of n&, the thermal distribution is: 

= 1 ( nb ) 
1 + n b VI + n b J 

As shown in Section 3 above, the asymptotic behaviour of the approach 
to equilibrium and the long time correlation functions are both determined 
by the leading eigenvalues of the matrix S. Fig. 1 shows how the correlation 
length Eq.(|T0"D changes, for n b = 3, with pumping rate and interaction time 
0. The axes correspond to the pumping rate N and 9 = gr\fN \ which is a 
scaled time parameter that is useful in revealing the phase structure. One 
can readily see the evidence for critical points at 9 ~ 1 and 9 ~ 4.6, which 
mark the transition from the thermal phase to the maser phase and from 
the thermal phase to the first critical phase, respectively. As anticipated the 
phase structure matches the one obtained by [|J using the continuous master 
equation. 

In the present work, we implement Eq. (|32|) directly by doing the numerical 
"experiment" of sending in one atom after another (i.e. multiplying p by 
S), and seeing how the photon distribution changes with k, the number of 
atoms that have passed through the cavity. The purpose of the numerical 
experiment was to measure how long it takes to get to equilibrium for different 
values of the physical parameters iV and n^. This could be important in 
physical experiments in which the results are interpreted in terms of the 
equilibrium photon distribution. In general we found that convergence was 
very rapid, with some interesting anomalies in the short time behaviour 

3 In order to calculate the eigenvalues we truncated the photon number at n = 200, so 
that S was a 200 x 200 matrix. 

4 In the subsequent analysis gr is kept fixed, which is relevant experimentally. The 
effect of varying gr will be examined in future work. 

5 Short "time" here actually refers to the first few atoms in the iteration. 
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(73) 



(74) 



In order to deal with finite dimensional matrices, we need to truncate 
the photon distribution at n — n max , say. Consistency requires that the 
probability of having n max photons in the cavity be small compared to the 
numerical accuracy of the calculation. This was checked by calculating the 
normalization of the probability distribution after each iteration. We found 
that the slight error in the normalization grew geometrically with the number 
of iterations, which was problematic for runs that contained thousands of 
atoms. We found however that this behaviour could be corrected by simply 
re-normalizing p(k) at each iteration. If this was done, the errors grew only 
linearly with k, and n max of about 200 was sufficiently large for our purposes. 

The purpose of the numerical experiment was to measure how long it takes 
to reach equilibrium for different values of the physical parameters gr, N and 
rib- This could be important in physical experiments in which the results 
are interpreted in terms of the equilibrium photon distribution. In general 
we found that convergence was very rapid. In order to have a quantitative 
measure of "how close" the system is to equilibrium, it is necessary to define a 
suitable measure on the space of photon distributions, which for the present 
purposes could be thought of as an n max dimensional vector space. We 
therefore define the distance between the photon distribution after k atoms 
and the equilibrium distribution by: 

Ap n (k) = \ Pn (k) - pT Mb \ (75) 

and take as a test for convergence the condition that all p n {k) must be within 
a certain range of the equilibrium value: 

max Ap n (k) < 0.005 (76) 

In order to have a point of comparison, we also checked convergence with a 
different measure, namely: 

AV + (k) = \V+(k) - V+(equil)\ (77) 

where V + (k) and V(equil) are the probability for an atom entering emerging 
from the cavity in the excited state for photon distributions during the transit 
given by p n {k) and p^, respectively, (cf. Eq. ([To|) .) We then used as a second 
test for convergence 

AV+(k) < 0.005 (78) 
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This test compares the probability that the kth atom will emerge in the 
excited state to the same probability at equilibrium. It therefore has direct 
physical relevance. Figs. 2a) and 2b) plot Ap(k) and AP + (k) as functions 
of k for N = 45, nj, = 3.0 and r = 1.0. It is interesting that the system first 
moves towards equilibrium (very rapidly for AP+(k)) and then moves away 
from equilibrium before it settles into its exponential approach to equilibrium. 
This feature appears to be fairly generic for large N. 

In order to do a systematic analysis of the convergence rate as a function 
of the physical parameters, we define k max as the number of atoms it takes 
for Ap(k) to get to some critical value, A crit , or less. For small enough A crit , 
k max should be large, in which case the convergence will be determined by 
the next to leading eigenvalue of S. Fig. 3 plots k max as obtained from the 
two tests as functions of 9 and N for rib = 3, with the condition A cr i t = 0.005. 
Clearly the resulting values of k max are large enough to be in the region of 
asymptotic convergence and the phase structure is the same as that predicted 
analytically using the eigenvalues of S. This result confirms the validity of 
our numerical method. 

We now use the numerical experiments to investigate a different but re- 
lated property of convergence. In particular, we look at how k max is affected 
by a change in rib arid N for fixed interaction time. As shown in Fig. 4, as rib 
gets large, the convergence is uniformly rapid for all N, whereas for low rib 
(i.e. low temperatures) there are critical values of N for which convergence 
slows suddenly. These are presumably the same transitions as in Fig. 1 but 
seen from a different view. In particular, lines of fixed rib and gr correspond 
to a section of Fig. 1 with 9 = gr^/N. Fig. 5 shows this structure for fixed 
rib and gr. As shown in Fig. 6, these "steps" in the convergence coincide (at 
least approximately) with values of the parameters in which the correlation 
length, as determined by the leading eigenvalue, increases. Fig. 6 plots the 
correlation lengths for n& = 1 for the first four leading eigenvalues. It is in- 
teresting that the "steps" correspond to points where the eigenvalues appear 
to cross. However, the eigenvalues are not degenerate 0, so the curves do 
not actually intersect. In Fig. 7 the equilibrium photon distributions are 
plotted for rib = 1 as a function of N. It shows clearly that the "crossing" 
of the eigenvalues is related to a discrete shift in the peak of the equilibrium 
photon distribution. The crossings, and the associated transitions in the 
convergence rate illustrated in Fig. 5, occur when the photon distribution is 
in the process of shifting, i.e. where there are two peaks. 
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5 Conclusions 



We have presented a systematic analysis of the discrete master equation de- 
scribing the approach to equilibrium of the micromaser in the presence of 
thermal dissipation. As expected, the long time behaviour is determined by 
the leading eigenvalues of the discrete transformation matrix. Interestingly, 
the eigenvalues of the matrix, evaluated numerically, become at some places 
nearly degenerate, and the phase structure of the micromaser occurs at or 
near where the eigenvalues come close to crossing. Our analytic results con- 
firm general features that emerge from the continuous master equation in || . 
We also have examined the approach to equilibrium of the micromaser both 
at short and long times using numerically methods. Our numerical results 
are consistent with the behaviour expected from the leading eigenvalues of 
both the discrete and continuous transformation matrices. However, our re- 
sults also show some interesting feature for short times; in particular, the 
system in general first approaches equilibrium relatively rapidly, then moves 
away from equilibrium, and then finally settles into its exponential approach 
to equilibrium. This behaviour appears to be fairly generic for large values 
of N, but more analysis is required to determine the source and relevance of 
these short time features. Moreover we have, in the numerical experiments, 
kept the transit time fixed. In future work we hope to examine how varying 
gr affects the above mentioned features. 
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Appendix: The Relationship Between the Continuous and Dis- 
crete Cases 

We expect differences in the dynamical behaviour in the discrete and 
continuous formalisms, but there should exist a limit in which the two for- 
malisms coincide. We look at the discrete formalism in the large flux limit. 
We take k = Rt to be large, which means that t ^> 1/T, or that the to- 
tal time over which the system is observed is much greater than the time 
between individual atoms. 

To take the large flux limit of the discrete master equation, we follow the 
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derivation of Eg. flB3|) . The discrete master equation has the form Eq.i 



We write, 







p k+l = Sp k 




(79) 


p k+l 


- p k 


= (s-i) P k 






(p k+1 - 


-p k ) 


= T[T-\S- 


l)T]T- 1 p k 




P k b +1 






l)^nm] (T V ) ma 


(80) 



which can be written in differential form as, 

dp=-(l-S)pdk (81) 

Since k = Rt we can write dk — Rdt and compare Eg . (j3ll) and Eq.(|8l|): 

- 7 L = -(l-S)R 

L = N(l-S); Rh = N (82) 

Using Eq. (p8|) and Eq.(^) we find, 

L c - N{M - 1) = N - N (83) 

1 + Lq/N 

We study the behaviour of the next to leading eigenvectors, since the corre- 
sponding eigenvalues control the behaviour of the approach to equilibrium. 
We have, 

(L c - N(M - l))pW = A<V 2) (84) 



and 



1 + L c /N 



SpF> 



K m p m 
K m p (2) 



(Lc-^^)PV = ^J-1]P (2) (85) 



Comparing Eq. (|8~4]) and Eq. (|85|) we have, 



N = N/k®; A< 2) = N[±r - 1] (86) 
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Solving this set of equations we have, 



_i N 

-R\n K < 2 > ~ ^\' 2 > = 7 A< 2 ' (87) 

where we have taken N large, or I/7 ^> T, which means that the typical 
photon decay time is much greater than the typical time between atoms. 
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Figures 




Figure 1: Variation of ccorrelation length with pumping rate and interaction 
time for fixed = 3 
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Figure 2: Variation of (a) Ap(k) and (b) AP + (k) with k 



22 




Figure 3: Variation of k max with 6 and N for fixed = 3 for the two tests 
described in the text. 
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Figure 4: Variation of k max with n& and N for fixed interaction time. 
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Figure 5 
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Figure 5: Variation of k max with N for fixed nj, and gr. 
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Figure 6 
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Figure 6: Variation of correlation lengths with N for rib = 1 for the first four 
leading eignevalues. 



26 




Figure 7: Plots of the equilibrium photon distributions for nj = 1 as a 
function of N 
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